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We calculate the ghost two-point function in Coulomb gauge QCD with a simple model vacuum 
gluon wavefunction using Monte Carlo integration. This approach extends the previous analytic 
studies of the ghost propagator with this ansatz, where a ladder-rainbow expansion was unavoid- 
able for calculating the path integral over gluon field configurations. The new approach allows us 
to study the possible critical behavior of the coupling constant, as well as the Coulomb potential 
derived from the ghost dressing function. We demonstrate that IR enhancement of the ghost corre- 
lator or Coulomb form factor fails to quantitatively reproduce confinement using Gaussian vacuum 
wavefunctional. 
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I. INTRODUCTION 

A combination of analytical calculations, based on 
Dyson-Schwingcr equations [1-13] and lattice gauge sim- 
ulations [14-36], has given new insights into the behav- 
ior of QCD Green's functions. In particular, it has been 
found that in the Landau gauge at low momentum the 
ghost propagator is enhanced while the gluon propaga- 
tor is suppressed. Dyson-Schwinger equations potentially 
admit solutions that are critical in the infrared (IR), i.e. 
the ghost propagator is divergent and the gluon propa- 
gator vanishes at zero momentum. On the other hand, 
the interpretation of lattice results is still somewhat con- 
troversial, since the IR region is sensitive to finite vol- 
ume effects and possible lattice artifacts in mapping be- 
tween the continuum and lattice definition of propaga- 
tors [2, 27-30]. One of the original motivations for such 
studies follows from the observation that for the phys- 
ical spectrum to consist only of color singlet states it 
is necessary that the ghost and gluon propagators are 
critical (in the sense defined above) [37, 38]. The ab- 
sence of colored states in the physical spectrum is often 
taken as a manifestation of confinement. The relation 
between the IR behavior of the ghost and gluon prop- 
agators and the expectation value of the color charge is 
tied to the realization of the residual gauge symmetry re- 
maining after imposing the Landau gauge condition [39] . 
The connection between remnant gauge symmetries and 
confinement, however, remains an unsettled issue, there- 
fore so does the relation between the IR behavior of the 
propagators and confinement. The relationship between 
the gluon and ghost propagators and confinement can be 
investigated in other gauges, and the Coulomb gauge can 
be particularly illuminating [40-45]. 

In the Coulomb gauge the time component of the 
vector potential becomes constrained by the trans- 
verse gluon field defined in the spatial directions alone. 



A a (x) satisfies, V • A a — for all color components, 
a = 1 • • • Nq — 1 , leading to an instantaneous poten- 
tial between color charges. This potential depends on 
the inverse of the Faddeev-Popov, or ghost, operator, 
M- 1 (A) = [V • B(A)}- 1 , with B(A) being the co- 
variant derivative in the adjoint representation. It was 
postulated by Gribov [46] and by Zwanziger [47] that 
gauge field configurations near the boundary of the field 
space domain, the Gribov horizon, dominate matrix el- 
ements and, since at the boundary the Faddeev-Popov 
operator vanishes, the instantaneous Coulomb potential 
is expected to be enhanced compared to the value at 
zero field [46-50]. This could signal confinement. Fur- 
thermore, since for a state containing a static quark- 
antiquark pair in the vacuum the Coulomb potential 
provides an upper limit on the total energy, Zwanziger 
concluded that a necessary condition for confinement is 
that the expectation value of the Coulomb potential in 
such state is also confining [40] . From the point of view 
that the energy spectrum is a direct probe of confine- 
ment, it seems relevant to investigate matrix elements 
of the inverse of the Faddeev-Popov operator. Analyt- 
ical calculations have been performed in, for example, 
Refs. [42, 44, 45, 51-55]. These typically start from 
an ansatz for the vacuum wave functional and various 
approximations are used to derive Dyson equations for 
correlations functions. Since the Coulomb energy in- 
volves fields at one time slice, only spatial correlations are 
needed. Through a systematic study of the IR behavior 
of the gluon-gluon correlation function and the Faddeev- 
Popov operator it was shown that within the particular 
set of approximations used to derive the Dyson equa- 
tions, all self consistent solutions are IR finite, but close 
to being critical. Most likely what this means is that the 
vacuum wave functionals used in these calculations do 
not yet account for all field configurations responsible for 
confinement. Another way of seeing this is through the 
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behavior of the spatial Wilson loops, for which such wave 
functionals fail to reproduce the area law behavior. If 
and when missing configurations are properly accounted 
for one would still face the question of reliability regard- 
ing the other approximations used in deriving the Dyson 
equations. These are typically based on the large- Nc ex- 
pansion and examination of the IR and ultraviolet (UV) 
behavior of higher order diagrams. To leading order this 
amounts to summing the rainbow-ladder diagrams. 

In this paper we confront the Dyson equations for 
the Coulomb gauge correlators with direct evaluation of 
the underlying matrix elements using Monte Carlo tech- 
niques for the path integral over the transverse gluon 
fields. The numerical techniques are close in spirit to 
those of lattice gauge theory, and are detailed in Sec- 
tion III. We begin by giving, in Section II, a short sum- 
mary of the Coulomb gauge and derivation of the Dyson 
equation. A summary and conclusions are given in Sec- 
tion IV. 



II. COULOMB GAUGE QCD 

In the Schrodinger representation the degrees of free- 
dom of the Coulomb gauge Yang-Mills theory arc: the 
transverse gluon fields, A a (x), which are the general- 
ized coordinates, and their conjugate momenta II a (x) = 
—i5/SA a (x.), equal to the negative of the transverse 
chromo-electric field [56]. These satisfy the canonical 
commutation relation, 

[IF' a (x),^> b (y)] - -i<W^'(V x )S(x-y), (1) 

where <5^ is the transverse projector (5^(V) = <5jj — 
ViVi/V 2 . The canonical Hamiltonian is a function of 
the generalized coordinates and momenta, and is given 

by 

H=^Jd* [j-^wj • rr(x) + B a2 (x)] + v, (2) 

where the chromo-magnetic field, B, is given by, 

B a (x) = V x A a (x) + |.f hc A b (x) x A c (x). (3) 

As usual, repeated indices are summed over. In Eq. (2), 
J — det(M(A)) represents the curvature of the Coulomb 
gauge field domain and is given by the Jacobian of the 
transformation from the ^4° = (Weyl) gauge - which 
has a flat field space - to the Coulomb gauge. Here, M 
is the Faddeev-Popov operator, 

M ab (x, y) = [-X7 2 J ab + gf abc A c ■ V x ] <5 3 (x - y). (4) 

The Coulomb potential, V, is obtained by using the equa- 
tions of motion to eliminate the longitudinal gauge field, 
and can be written 

V = \ J d^yj- V a (x)^ afc (x, y; A)p b (y), (5) 



where, in the absence of quarks, the color charge density 
is given by 



p a ( x ) = / a6c n"(x) • A c (x), 



(6) 



and the Coulomb kernel , K (A) is 

K(A) = gM-\A){-V 2 )ghr\A). (7) 



In the abelian limit this kernel reduces to, 

g 2 6 ab 



47r|x-y| 



(8) 



the familiar expression for the Coulomb potential be- 
tween charges located at points x and y. Denoting 
the vacuum wave functional by ^f[A] = (A\^), the vac- 
uum expectation value, (vev) of an operator 0[A] in the 
Coulomb gauge is given by, 



(0) = 



(*\Q\V) 



(9) 



where 



<tf|0|tf) = / VAJ[A]0[A]\^[A]\ 2 , (10) 
Ja 

and the integral is restricted to the fundamental modular 
region (FMR) A e SI which is inside the Gribov region 
O. The FMR is defined as the set of gauge fields A a (x) 
corresponding to the absolute minima of the functionals 
I[g] = J e?x(A a9 (x)) 2 minimized with respect to time- 
independent gauge transformations g — g(x), while the 
Gribov region £1 also includes local minima of /. It has 
been argued by Zwanziger [57] that the bulk of the inte- 
gral measure is concentrated on the common boundary 
of FMR and the Gribov region and in the Monte Carlo 
simulations presented here only the restriction to ft will 
be implemented. The vev of the inverse of the Faddeev- 
Popov operator, which in the Coulomb gauge plays the 
dual role of the ghost propagator and the running cou- 
pling, is given by 



d(k) 



1 5 ab f 

N 2 -l J 



dxe 4kx (*| ff M^ 1 ' a6 (x,0)|*), (11) 



where d(k) is referred to as the ghost dressing function; 
at tree-level, d(k) = 1. If the expectation value of the 
Coulomb kernel is approximated by the square of the 
vev of the ghost propagator then the momentum space 
Coulomb potential between a color-singlet static quark- 
antiquark pair becomes V(k) = -Cpd 2 (k)/k 2 [41, 42, 
44]. In general, however, one expects the two vevs to be 
different and this difference can be accommodated via an 
additional form factor and results in the potential of the 
form V(k) = -C F d 2 (k)f(k)/k 2 [51, 58]. It is clear that 
if the ghost becomes IR enhanced, d(k) » 1 as k — > 0, 
the Coulomb interactions between color charges becomes 
stronger as the separation between charges increases. To 
obtain a linearly rising potential, however, it would be 
necessary for the product d 2 (k)f(k) to be critical with 
d 2 (k)f(k) -> k~ 2 ask->0. 
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A. Dyson equations 

The set of coupled Dyson equations for the ghost dress- 
ing function d(k), the Coulomb dressing function j(k) 
and the gap equation, which determines the gluon-gluon 
correlation function, were derived and extensively stud- 
ied in Refs. [45, 51-55]. Here we only summarize the 
main features of the ghost and gluon correlation func- 
tions. In these studies the vacuum wave functional was 
parametrized as a gaussian 



%?[A] = exp 



d 3 k 

(2^)3 



,(fc)A a (k)A a (-k) , (12) 



with uo(k) being a parameter. It was shown in Refs. [52- 
55] that, to leading order in the loop expansion, the ef- 
fect of the curvature J could be absorbed by a redef- 
inition of uj with the gap equation correlating the low- 
mometum behavior of u> and the curvature. In the sub- 
sequent derivations of the Dyson equations we thus set 
J = 1 The vacuum wave functional can be optimized 
by minimizing the vacuum energy density with respect 
to oj(k). This leads to a gap equation which after renor- 
malization depends on the renormalized coupling g r {^) 



and the boundary condition uj{k — ► 0) 



As long as 



m g is finite one finds that the solution of gap equation is 
qualitatively insensitive to g r (fj) and can be well describe 

by, 



u>(k) 



m g if k <m g 



otherwise. 



(13) 



It should be noted that m g is a mass parameter intro- 



mit) 




FIG. 1: Diagrammatic representation of the expansion of the 
functional integral for the ghost propagator c.f. Eq. (14) 

duced by the ansatz wave functional and should not be 
confused with the QCD scale introduced by renormaliza- 
tion. The latter appears in the renormalized Dyson equa- 
tion for the ghost dressing function which, as mentioned 
earlier, can be identified with the running coupling. In 
principle, m g = m g (g r (u,), fi) should be renormalization 
point invariant and just like g r (fJ-) determined by a phys- 
ical observable e.g. the string tension. Within the set of 
truncations build in the derivation of the Dyson series, 
most likely the renormalization group invariance of m g 
can not be proven and we shall consider m g as a free pa- 
rameter. Given u>(k) the Dyson series for the ghost dress- 
ing function can be sum up and represented as a single 



integral equation within the rainbow-ladder approxima- 
tion, illustrated in Fig. 1. All omitted diagrams have 
at least one vertex loop correction (e.g. last diagram in 
Fig. 1), which were shown to be generally smaller than 
the self-energy loops [51]. The diagrams shown in Fig. 1 
represent functional integrals over | , FL4]| 2 of polynomi- 
als of the A field originating from the expansion of the 
inverse Faddeev-Popov operator 
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(*l*5 
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-V 2 



gMA}\ 2 
1 .' . / VAgA x 



— gA x — 
V 2 V 2 



(14) 



where x refers to the color space. Neglecting the re- 
striction to the Gribov region enables one to perform 
the functional integrals analytically, and neglecting con- 
tractions that corresponds to vertex corrections makes it 
possible to re-sum the series, resulting in, 
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1 



d(k) 5 (A) 



l-(k-q) 2 



(15) 

The dependence of the bare coupling, g = g(A), and 
the loop integral on the UV cut-off has been shown ex- 
plicitly. Instead of using the bare coupling and the UV 
cutoff as the renormalization point, the equation can be 
renormalized at a finite momentum scale through sub- 
traction, which also defines the renormalized coupling as 
g r (fi) = d(k = u) 



1 



1 



d(k) 



N, 



c 



d{u) 
<iq 



l-(k-q) 2 ri(|k-q|) 
(k - q) 2 2co(q) 



M) 
(16) 



As discussed above, the mass scale is brought in through 
the function uj, and in the case discussed here, it is given 
by m g = w(0). Thus from now on we will use the nota- 
tion k = k/m g to denote dimensionless momenta. The 
solution of the Dyson equation for the ghost propagator 
depends on one more parameter, the value of d(p) at a 
single point, i.e. at p, = u,jra g = 1. In Fig. 2 we plot 
the numerical solutions of Eq. (16), as a function of mo- 
mentum in units of m gi for three choices of d(k =1). As 
d(l) is increased the solutions become more IR enhanced 
until at, approximately, d(l) ~ 3.41 the solution becomes 
critical [51]. Above this critical point the Dyson equation 
has no solutions, i.e. develops a Landau point at physi- 
cal, k > momentum. This is a sign that the functional 
integration in Eq. (14) has crossed the Gribov horizon. 
The mass scale dependence of the ghost propagator can 
be best understood by using an angular approximation, 
to the integral in Eq. (16), 



6{k- q)k + 6(q-k)q, 



(17) 
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FIG. 2: Comparison between the numerical solutions of the 
Dyson equation for the ghost propagator and the approximate 
analytical solutions of Eqs. (18), (19). The solid (dashed) lines 
represent the low (high) momentum behavior, respectively. 



which enables one to transform the integral equation to 
a first order differential equation that can be solved an- 
alytically and further well approximated by [51], 



d(jfe) 
d(k) 



d(jl) 



l + /? L dV7(£) (fc-/i) 

<m 

l + /fedV7(#io g (!) 



T , k<l, (18) 

y,k>l, (19) 



where 7 = 1/2 and /3 L = (5/6)(7V c /3)/7r 2 and (3 H = 
(N c /3)/ir 2 i.e. L ~ /3 H - 0.1 for 7V C = 3. It clearly 
follows that the ghost propagator is independent of the 
renormalization scale, /x and depends on a single value of 
ci(/i) at an arbitrarily chosen renormalization point. Fur- 
thermore, from Eq. (18) it follows that a solution exists, 
i.e there is no Landau pole, as long as d(p) < 1/(/3l/x) 7 . 

As discussed above, the approximations leading to 
Eq. (16) include eliminating all vertex corrections and 
neglecting the restriction on the functional integral to be 
contained within the Gribov horizon. In the following 
we present results from a Monte Carlo simulation of the 
ghost propagator that does not have these limitations. 



III. MONTE CARLO CALCULATION 

The evaluation of the functional integral in Eq. (10) 
is usually performed analytically by expanding the op- 
erator O in a power series over the gauge field A and 
truncating at some order (c.f. Eq. (14) ). Here we avoid 
these approximations by evaluating the functional inte- 
gral by Monte Carlo integration using the model wave- 



function (12) with the approximate solution (13) for 
w(fc). 

The gluon configurations are generated on a = 3 di- 
mensional momentum space grid. The gluon fields A^(k) 
are N4x(Nq — 1) complex numbers per lattice site, where 
Af(—k) = Af*(k). The momentum is discretized on the 
lattice as 



ki 



ie {1,2,3}, 



(20) 
(21) 



where a denotes the lattice spacing. The gauge fields 
must satisfy the position space Coulomb gauge condition 



Y,A a i {x)-A a i {x-ai)=0, 



(22) 



i=l 

which translates in the momentum space to 
3 



^(1 - cos(afc) + i s\n(ak))A a t (ak) = 0. (23) 



From now on we will use the notation A{k) = A(ak)/a 2 , 
etc. in reference to dimensionless quantities scaled with 
the lattice spacing. The coupling is incorporated by gen- 
erating gAf(k) rather than Af(k), which requires sub- 
stituting u>(k) with u){k)/g 2 in the model wavefunction. 
The gluon fields are generated with the distribution 



MA]f 



exp 



1 



EE E Mmi-k)"f 

m i=i 0=1 y 



(24) 

This is accomplished by independently generating two of 
the vector components, Af, with a heatbath, then con- 
structing the third component such that the momentum 
space Coulomb gauge condition, Eq. (23), is satisfied. 

The calculation of the Jacobian is akin to the calcula- 
tion of the quark determinant in lattice QCD and in the 
present work it is set to one. The Jacobian was included 
in Ref. [52] in a certain truncation scheme. There it was 
found to lessen the dependence of the Coulomb potential 
to the choice of coupling. 

As a first test we evaluate the gluon propagator, 



g 2 G(k) 



1 1 

NlN d -lN 2 



N d -lN 2 c -l 

TjfEE A«Ck)AU-k)). 



i—1 a—1 



(25) 

The value of G(k) is analytically known to be G(k) — 
l/2w(fc). The numerical result, shown in Fig. 3 does in- 
deed agree with the analytical one, where the numerical 
statistics are improved by taking the Z 3 average, that is, 
averaging over the three equivalent directions in momen- 
tum space. Since k = m g k The physical propagator in 
units of m g is given by 



m g G{k) — rhgG(k). 



(26) 
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FIG. 3: (Color Online) The gluon propagator calculated on 
1000 gauge field configurations with various parameters. The 
analytic results are plotted as a continuous line, showing a 
good agreement with simulations. 

We now proceed to computing the ghost dressing func- 
tion. The ghost propagator is expressed as the expec- 



tation value of the inverse of the Faddeev-Popov (FP) 
operator, Eq. (11). The discrete form of the FP operator 
was derived in Ref. [47] 

Nd t \ 
M ab (x, y)4>\y) = ]T S ab (<j> b (y + i) + $ b (y - i) - 2^{y)) 

i=l 

-f abc (4> b (y + t)Miv) " 4>"(y - t)A c S - *)) , 

which is real and symmetric. Note that the region of 
integration in Eq. (10) is the Gribov region where M 
is positive definite. Thus any gauge field configuration 
that produces a Faddeev-Popov operator with negative 
eigenvalues must be discarded. 

With periodic boundary conditions imposed on the lat- 
tice, M ab (x,y) has - 1 trivial zero modes, making it 
formally non-invertible. This problem is avoided by fol- 
lowing Ref. [22] and solving 

jf dyA/ Qb (x, y)0 b (y) = S ab (<5(x) - I) . (27) 

The position-color vectors are then Fourier transformed 
to momentum space and the inverse of M ab (a;,0) recov- 
ered. 



dxe 



-ikx 



t> a (x)) = 



£k-x 



((M^r^o))- 



V 



dxdye 



-?k-x 



((M- 1 r(x,y)> = ^-5(k) / 
9 Jv 



.x^M, 

9 

(28) 



where D(k) = d(k)/k 2 . The average is taken over gauge 
field configurations, and finally the ghost propagator is 
Z$ averaged. In the free case, M — » —V 2 and the prop- 
agator would be 

t)(k) = 9 —, = A, (29) 

4E 4 sin 2 (fc,/2) fc 2 ' 

which defines an appropriate momentum variable. The 
same philosophy is used in conventional lattice QCD 
studies of gluon propagator [16, 59, 60]. 

With the model wavefunction (12), the coupling g is a 
free parameter. The larger g is chosen to be, the broader 
the Gaussian. This increases the fluctuations of the gauge 
fields and M develops smaller eigenvalues, resulting in 
the infrared enhancement of (M~ 1 (fc)). With increas- 
ing the value of g, the FP operator becomes likely to 
develop negative eigenvalues. While this means that a 
(possibly large) proportion of the generated gauge fields 
must be rejected, it is necessary for the entire domain of 
the functional integration to be sampled. The number of 
rejected configurations grows rapidly when the value of g 
approaches certain critical value, which depends on the 
value of m g used in the model for uj(k) of Eq. (13). This 



is easily understandable, as larger value of m g means the 
gluon wavefunction is infrared enhanced in a larger inter- 
val of momenta, yielding narrower Gaussian width over 
that interval. 

Each generated gauge configuration used in calculat- 
ing the ghost dressing function is checked to lay in the 
Gribov region by calculating several eigenvalues of the 
FP operator to ensure their positivity. If the latter con- 
strain is not imposed, the resulting ghost propagators are 
dominated by numerical fluctuations (resemble random 
noise) in the region where the generated gauge configura- 
tions have large fraction laying outside of Gribov region. 
For example, for rh g — 1.5, the fraction of rejected con- 
figuration ranges from nearly 0% for g < 1 to 100% for 
g > 1.1 with sharp increase above g = 1. For rh g = 5 this 
"critical" value of g increases to about 1.6 . In our calcu- 
lations we restrict to the region of g, where the fraction 
of rejected configurations does not exceed 20% to main- 
tain moderate computational time. The resulting ghost 
dressing function is shown in Fig. 4 for a calculation with 
1000 gluon configurations on a 40 3 lattice with rh g = 1.25 
and g — 0.7. 

In order to relate the calculated ghost dressing function 
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to the physical region several issues should be resolved 
that would allow to draw a correspondence. Here we 
review the most relevant ones. 



A. Lattice Artifacts 

Discretization of space introduces several artifacts, 
that should be accounted for. These are errors introduced 
by finite lattice volume, finite lattice spacing, which also 
induces broken spatial rotational symmetry. 



1. Finite Lattice Spacing 

It is argued in the Refs. [16, 59, 60] that using the 
redefined lattice momentum variable of Eq. (29) allows 
one to avoid the leading-order discretization errors aris- 
ing from the ultra-violet cutoff in momentum introduced 
by the finite lattice spacing. Still, the errors from reduc- 
ing the spatial rotational symmetry 0(3) group down to 
discrete Z(3) are unaccounted for. These manifest them- 
selves as a large spread in the ghost propagator. This 
spread occurs in a characteristic pattern, as can be seen 
in Fig. 4, which becomes more prominent with increased 
lattice volume. These patterns can be easily understood 
by considering a selection of subsets of the points plotted 
by using certain criteria imposed on the momentum vari- 
able. The first subset considered has the constraint that 
all three components of the momentum are equal to each 
other (laying on the diagonal direction of the lattice). 
This selection of the data forms a smooth line through 
the upper part of the plot. A subset including points 
with two of the momentum components equal to each 
other and the third one set to zero (along the diagonal 
direction of the cube's side) forms another smooth curve, 
this one going through the middle of the plot. Finally, 
the subset with only one non-zero component of momen- 
tum (along the side of the cube) forms a line passing 
through the lowest part of the plot. These subsets are 
shown in Fig. 5a. Furthermore, if the constraints de- 
scribed above are allowed to be violated by a few units 
of minimum lattice momentum, the rest of the points in 
the plot start to fall into these subgroups, as shown in 
Fig. 5b. Thus, for the further analysis of our data we will 
use only a subset of points with momentum components 
not differing from each other by more than one unit of 
minimum lattice momentum. This is the "cylinder cut" 
introduced in Refs. [16, 59], which allows us to select the 
points least affected by errors introduced by the broken 
rotational symmetry and leaves a sufficient number of 
points for statistical analysis. 



2. Finite Volume Effects 

While a consistent treatment of the finite volume ef- 
fects requires extensive investigation into discretization 
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FIG. 4: (Color Online) Ghost dressing function d(k/m g ) ver- 
sus k/m.g for rh g = 1.25 and g — 0.7 on a lattice with 40 3 
volume. Here the crosses denote the full data set and open tri- 
angles denote the subset of points with equal momenta com- 
ponents. 



of the theory on the lattice, here we simply investigate 
this dependence by comparing benchmark calculations 
on lattices with different volumes. A set of calculations 
with four different lattice volumes is shown in Fig. 6, 
which shows that there are very small variations only 
in the low momenta region for lattice volumes from 20 3 
to 40 3 . Thus we choose to use lattice volume of 20 3 for 
the further calculations, which allows for both reasonable 
computational time and small errors. 



B. Renormalization 

The introduction of a finite momentum grid provides 
a sharp cutoff for the regularization of the ultra-violet 
divergences, i.e. it is equivalent to the role of A in 
Eq. (15). In order to identify the ghost dressing func- 
tion with the running coupling, for each lattice spacing 
it should be possible to choose the lattice coupling, g in 
Eq. (24), so that the results of simulations are indepen- 
dent of the lattice spacing. In the simulation, explicit 
dependence on the lattice spacing enters through depen- 
dence on rhg = am g , e.g. when the lattice ghost dressing 
function is plotted against k = k/m g = k/rh g the result 
should be independent of rh g and depend only on the 
value of the renormalized coupling. In practice, we pro- 
duce a series of simulations with different values of g in 
the range between g = 0.5 — 1.5 and rh g in the range of the 
accessible lattice momenta 2it/NL at < rh g < V37T. We 
then compare with the scaling predicted by the solutions 
of the Dyson equations given in Eqs. (18) and (19), in the 
low and high momentum region, respectively. In the high 
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FIG. 5: (Color Online) Ghost dressing function for d(k/m g ) 
versus k/m g for rh g = 1.25 and g = 0.7 on a lattice with 40 3 
volume. Here the crosses denote the data points and the three 
lines connect the subsets of the points laying within a) and 7 
b) minimum lattice momenta correspondingly of the edge, the 
side diagonal and the diagonal directions of the momentum 
lattice cube. 



momentum regime, k/m g — k/rh g is kept large by run- 



> 



2n/N i 



Lat • 



ning simulations with small rh g i.e. with 
In this regime the constituent gluon mass is close to the 
minimum accessible momentum scale on the lattice and 
the gluon propagator is close to asymptotic while the 
non-perturbative effects are only present for a few, low- 
est momentum points. This regime should be described 
byEq. (19), 



d(k) 



d(ft) 



/3^/7'(/i)log(|) 



(30) 



V=12", NCont=1000 
V=2(f, NConkWOO 
V=3(f, NConf=1000 
V=40 3 , NCont=300 




FIG. 6: (Color Online) Ghost dressing function d(k/m g ) ver- 
sus k/m g for rh g = 1.25 and g = 0.7. We show the influ- 
ence of varying the lattice volume V on the calculated data 
points with momentum components not differing by more 
than two units lattice momentum, as described in the text. 
Here NConf denotes the number of the sampled gluon con- 
figurations. 



where (3' H and 7' will be treated as fit parameters. For 
Nd = 20 each data set has 30 momentum points after 
the imposed "diagonal" cut described above. We choose 
6 data sets with rh g £ [0.1,1] and g £ [0.3,0.75], where 
the 25 highest momentum points can be considered to be 
in the asymptotic region. For each value of the coupling, 
g, the value of d(ft) is fixed by the data itself with ft 
set equal to the momentum cut-off, ft = y/3n /rh g . The 
formula in Eq. (19) is fitted to all 150 data points by 
varying j3' H and 7'. The resulting remarkably good fits 
are shown in Fig. 7 with the best-fit value of (3' H — 0.86(2) 
and 7' = 0.5(2). The data deviate from the perturbative 
form at intermediate momenta, which is to be expected. 

On the other hand, in simulations with large rh g i.e 
for fh g < VSir, Eq. (18) should apply. Then the con- 
stituent gluon mass is close to the largest accessible mo- 



mentum scale on the lattice. This regime is dominated 
by non-propagating gluons induced by non-perturbative 
dressing. Here we expect, 



d(k) = - 



d(ft) 



1 + {3> L dvy (ft) (h - ft 



(31) 



In this region we select a total of 8 data sets composed 
of rh g S [3,5] and g £ [1,1.4], where the 15 lowest 
momentum points can be considered to be in the non- 
perturbative region. Here, d(ft) is obtained from each 
data set itself, at ft chosen, to avoid finite- volume effects, 
to be the second lowest momentum point. In the low mo- 
mentum range a total of 120 data points was fitted vary- 
ing (3' L while keeping 7' = 1/2 which was previously de- 
termined from the high momentum fit. A sample of data 
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FIG. 7: (Color Online) Simultaneous fits to the ghost dressing 
function d(k/m g ) versus k/m g for m g > 2n/NLat- 
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ancies in the higher momentum region are expected as a 
consequence of deviations from purely non-perturbative 
behavior set by the asymptotic tail of the gluon propa- 
gator. 



IV. CONCLUSIONS 

We have computed the ghost correlation function by 
direct Monte Carlo simulation of the functional inte- 
gral with a model gaussian wave functional and com- 
pared it with the solution of the corresponding Dyson 
equation. We have found that the scaling behavior of 
the solution of the Dyson equation is reproduced in the 
simulation. This confirms that the corrections to the 
rainbow-ladder approximation are both IR and UV fi- 
nite, and do not change the scaling properties. The (3 
function obtained from simulations is, however, an order 
of magnitude larger than the one from the Dyson equa- 
tion. This is to be expected, since the Dyson equation 
does not properly take into account the boundary of the 
field space integral, and thus is expected to overestimate 
the magnitude of the allowed field values and thus of the 
critical coupling. The Monte Carlo simulation still needs 
to have the Faddeev-Popov Jacobian implemented, but 
that is not expected to qualitatively change the results. 

In our simulations we have found that positivity of 
the Faddeev-Popov operator is not sufficient to produce 
critical behavior. This needs to be investigated further, 
in particular on larger volumes; nevertheless, since the 
simple gaussian vacuum wave functional does not probe 
topological configurations {e.g. of magnetic disorder) it is 
not too surprising that the IR enhancement of the ghost 
correlator or Coulomb form factor fails to quantitatively 
reproduce confinement. For this purpose a wave func- 
tional of the type proposed in Ref. [61] should be tried. 



FIG. 8: (Color Online) Simultaneous fits to the ghost dressing 
function d(k/m g ) versus k/m g for m g < v37T. 



points with the corresponding fits are shown in Fig. 8 for 
the best fit value of j3' L = 0.81(2). Again, the discrep- 
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